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ABSTRACT 

A collisional A-body simulation using NBODY5 on a single CRAY YMP processor is 
followed well into the post-collapse regime. This is presently one of the largest particle 
numbers of all such models published, but some data for an even larger A produced by 
using special-purpose computers have recently been presented. In contrast to previous 
ensemble-averaged A-body simulations the noise here is low enough to just compare 
this one single run with the expectations from statistical models based on the Fokker- 
Planck approximation. Agreement is as good as could be expected for the case of 
the evolution of the Lagrangian radii, radial and tangential velocity dispersions and 
various core quantities. We briefly discuss approximate models to understand number 
and energy of escapers and the question of gravothermal core oscillations; although the 
system exhibits post-collapse oscillations they turn out to be directly binary driven 
and we cannot prove the existence of a gravothermal expansion at this particle number. 
Finally in a detailed examination of the wandering of the density centre we find in 
contrast to some previous studies a clear long-time period of the order of about 14 
half-mass crossing times. 
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1 INTRODUCTION 

One of the grand challenges of theoretical astrophysics is to 
understand the dynamics of globular star clusters. Apart 
from the intrinsic interest of understanding the behaviour 
of large A-body systems, the importance of the problem 
stems from its relation with current observational research, 
and from the aim to study the behaviour of even larger sys- 
tems such as galactic nuclei. Without understanding globu- 
lar clusters by means of special tools for investigating their 
evolution, the possibility of unravelling the mysterious be- 
haviour of galactic nuclei seems very remote. Unfortunately, 
the direct simulation of such rich stellar systems as globular 
clusters with star-by-star modelling is not yet possible. The 
gap between the largest useful A^-body models with N being 
of the order of 10 4 particles and the median globular star 
cluster (N ~ 5 x 10 ) can only be bridged at present by use 
of theory. There are two main classes of theory: (i) Fokker- 
Planck models (henceforth FP model), which are based on 
the Boltzmann equation of the kinetic theory of gases (Cohn, 
Hut & Wise 1989, Murphy, Cohn & Hut 1990), and (ii) 
isotropic (Lynden-Bell & Eggleton 1980, Heggie 1984, Bet- 
twieser & Sugimoto 1984) and anisotropic gaseous models 
(henceforth AG models, Louis & Spurzem 1991, Spurzem 
1992, 1994), which can be thought of as a set of moment 
equations of the Fokker-Planck equation. 



These simplified models are the only detailed models 
which are directly applicable to large systems such as glob- 
ular clusters. But their simplicity stems from many approxi- 
mations and assumptions which are required in their formu- 
lation. Examples are the assumption of spherical symmetry, 
which is inconsistent with the asymmetry of the galactic 
tidal field, or statistical estimates of the cross sections for 
the formation of close binaries by three-body or dissipative 
(tidal) two-body encounters, and for their subsequent grav- 
itational interactions with field stars. Such processes play a 
dominant role in reversing core collapse of globular clusters, 
which otherwise would inevitably lead to a singular den- 
sity profile with infinite density at the centre (Lynden-Bell 
& Wood 1968, Bettwieser & Sugimoto 1984, Elson, Hut & 
Inagaki 1987, Hut 1993). 

Recent detailed comparisons have established a good 
relationship between such simplified models and direct N- 
body simulations, but for the even more idealized case of a 
single mass system (identical individual mass in contrast to 
real star clusters) and rather low particle numbers (Giersz 
& Spurzem 1994, henceforth GS, for N = 250, 500, 1000 
and 2000; this paper already contains some information on 
the present N = 10000 direct model, and also Giersz & 
Heggie 1994a, b, henceforth GHI, GHII). Only in one case 
has a two-mass system been compared with direct A-body 
simulations (Spurzem & Takahashi 1995). Unfortunately 
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larger N post-collapse star clusters cannnot be understood 
easily by extrapolating the results for low N. Whereas two- 
body relaxation scales with N/logjN (with 7 w 0.11), or 
in other words, all TV-body results before core bounce (the 
point where binary effects reverse core collapse) match ex- 
actly if one scales the time with this factor (GHI, GHII), the 
formation of three-body and tidal two-body binaries intro- 
duces other timescales into the problem, which do not scale 
with the same factor of N. 

Therefore it remains very important to get more data 
from direct iV-body simulations with N as high as allowed 
by the present-day computational technology and software, 
in order to strengthen confidence in the statistical models 
based on the Fokker-Planck and other approximations. To 
follow for example an N = 10000 star cluster to core bounce 
one needs to calculate some 900 initial half-mass crossing 
times. To do this accurately a high-order integration scheme 
with high intrinsic accuracy and no artificial softening of 
the potential (as is very often used in so-called tree codes, 
Barnes & Hut 1986, for collisionless stellar dynamics, typi- 
cally models lasting some 10 or 20 initial crossing times) is 
necessary. 

The required high accuracy integration scheme NBODY5( 
(Aarseth 1985) includes a fourth order Adams-Bashforth- 
Moulton integrator based on force interpolation polynomi- 
als at four past time points. It includes as well a two- and 
more body regularization technique (for two-body see Kus- 
taanheimo & Stiefel 1965) to deal efficiently with close en- 
counters and a two-level timestep scheme (Ahmad-Cohen 
neighbour scheme, Ahmad & Cohen 1973), which updates 
the irregular force due to neighbour particles in much shorter 
time intervals than the regular, less fluctuating force exerted 
by the more distant particles. 

NBODY5 has recently been improved to NBODY6, 
which uses a Hermite interpolation for the force and its time 
derivative (based on only two time points) and a hierarchical 
block time step scheme (variation of time steps only allowed 
by a factor of 2 to keep as many time steps commensurate as 
possible, Makino & Aarseth 1992). The latter is suitable for 
implementation on parallel machines, and actually has been 
implemented in a slightly advanced version NBODY6H — h on 
parallel computers (Spurzem 1995). 

On the other hand special-purpose computers have been 
built to calculate forces for gravitational iV-body systems 
effectively (Sugimoto et al. 1990). There are also so-called 
field programmable arrays, hardware devices which can be 
programmed in a similar way by hardware switches (Brown 
et al. 1992), which may be used for direct ./V-body simula- 
tions in the near future (D. Merritt, pers. communication). 
Such computers became increasingly available in the past 
(for example the HARP-2 and HARP-3 computers opera- 
tional in Kiel and Cambridge, respectively) and in autumn 
1995 a Teraflop machine consisting of many such HARP- 
boards has been presented during the IAU symposium 174 
in Tokyo (Hut & Makino 1996), which would be able to fol- 
low core collapse of a reasonably realistic globular cluster. 
By using existing pieces of this machine it already became 
possible to integrate systems larger than N = 10 4 past core 
bounce (N = 3.2 • 10 4 , Makino 1996) and observe the onset 
of gravothermal oscillations. 

In this paper we report on one of the largest published 
collisional iV-body simulations so far, which was performed 



during the last two years on a single CPU of a CRAY YMP 
computer by using NBODY5. Although it was not possi- 
ble to detect gravothermal oscillations without any doubt in 
that run, we demonstrate that 10 4 is a large enough particle 
number to exhibit a predictable evolution of the system as a 
whole. Predictable means here, that there is already reason- 
able agreement of this one simulation with the expectations 
from statistical models based on the Fokker-Planck approx- 
imation. This increases confidence that it will be possible 
to detect gravothermal oscillations in even larger simula- 
tions, which are also predicted by the FP models. On the 
contrary, for lower N, to check compatibility with theoreti- 
cal models it is necessary to average a sufficient number of 
statistically independent iV-body models (GHI, GHII, GS). 
Single systems of N = 1000, for example, behave practi- 
cally unpredictably because of the inherent chaotic nature 
of the orbits in an encounter-dominated star cluster. For 
large N, however, such chaotic behaviour of the individual 
orbits is balanced by the decrease of the fluctuations, be- 
cause the relative energies of individual stellar encounters 
decrease compared to the total energy of the system, as N 
increases. 

After some technical remarks for the model simulation 
(Sect. 2) we report in the following on the time evolu- 
tion of some relevant observables in the iV-body simula- 
tion of an isolated, single point mass, N — 10 4 star cluster. 
These are mainly Lagrangian radii, velocity dispersions and 
anisotropy, core parameters (Sect. 3.1), and the results are 
compared with anisotropic gaseous models. Escaper proper- 
ties, the individual evolution of hard binaries and their role 
for the generation of post-collapse oscillations are briefly ex- 
amined (Sect. 3.2), and the movement of the density centre 
is extracted from the data and discussed in some detail (Sect. 
3.3). 

2 THE SIMULATION 

A 10 4 -body realization (equal masses) of Plummer's model 
was integrated by using the standard NBODY5 scheme on 
one processor of a CRAY YMP machine for 2760 iV-body 
time units. In the following we use standard iV-body time 
units (Heggie & Mathieu 1986), in which the total energy 
of Plummer's sphere is E = —0.25, G = 1, and M = 1 
(G, M, gravitational constant, total mass of the system, re- 
spectively). In such a system the individual body's mass 
is m = 1/N for an equal mass cluster, where N is the to- 
tal particle number. The initial half-mass crossing time is 
icro = 2V2 w 2.828. Hence we have integrated the sys- 
tem over a little less than 1000 i cr o, which took about 1500 
CPU hours on the CRAY, equivalent to approximately two 
months continuous running time. The accumulated error of 
the entire simulation in total energy was A_E = 8.08 • 10~ 3 . 
The accuracy of the integration was checked by the crite- 
rion that the relative change in energy within a given time 
interval (0.5 t„ ) should be less than 10~ 5 . Since NBODY5 
contains the Ahmad-Cohen neighbour scheme each particle 
has two timesteps, a short one, after which the irregular 
force from the neighbours (typically less than y/N particles) 
is updated, and another regular timestep (which is usually 
about ten times larger than the irregular timestep), after 
which the total force exerted by all particles is updated, 
and changes in the neighbour list are maintained. Such a 
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neighbour scheme saves a considerable amount of CPU time 
as compared to a scheme which at any step always calculates 
the full force, at least as one considers a non-parallel general 
purpose computer. Both the regular and irregular timestep 
are fully adaptive according to the indidvidual force vari- 
ations on the particle (for details see Aarseth 1985). Our 
model performed in total approximately 13.5 billion irregu- 
lar steps, 1 billion regular steps, and 500 million regularized 
steps for closely bound pairs, whose internal motion is in- 
tegrated separately in regularized coordinates. This means 
on average each particle was moved 1.35 million times (its 
regular force updated ten thousand times). A more quanti- 
tative information can be found on Figs, la, b, which show 
the CPU time used and the number of regular, irregular, 
and regularized (KS) steps per iV-body time unit as a func- 
tion of time (Fig. la) and the same data per time unit per 
particle (Fig. lb). The fluctuations after t = 2300 are due 
to the onset of the activity of extremely hard binaries in 
the core (core bounce is at t — 2380, see below). These bi- 
naries suffer from strong superelastic scatterings with field 
stars in the core, so the timestep distribution of the field 
stars is affected itself. It is remarkable, however, that from 
the beginning until the core bounce, when there is a strong 
core-halo structure, the CPU time rates do not increase by 
more than half an order of magnitude, which is a proof that 
the two-level neighbour scheme of NBODY5 is indeed very 
well tailored to simulate centrally concentrated systems. 



Figure la. CPU time in seconds, number of irregular, regular, 
and KS steps per iV-body time unit, as a function of time. 



3 RESULTS 

3.1 Spatial Evolution 

In the following we compare our iV-body results with those 
of a standard AG model (Louis & Spurzem 1991, Spurzem 
1992, Spurzem 1994, Spurzem 1996), using for the latter 
those parameters which gave best agreement between all dif- 
ferent kinds of star cluster models in previous work, namely 
A = 0.4977, Aa = 0.1 and C b = 90. Here A is the numer- 
ical parameter of order unity scaling the heat flux, related 
to the standard C constant in isotropic gaseous models (see 
e.g. Heggie & Stephenson 1988) by 

x = l^c. 

Our A corresponds to the standard C = 0.104 of Heggie 
& Stephenson (1988), which gives best agreement between 
FP and AG models. Aa is a numerical factor scaling the 
anisotropy decay time Ta, defined such that the decay of 
anisotropy by two-body relaxation is 

(SA/8tU = -j^r , 

where Ta is the self-consistent anisotropy decay time for a 
Larson type anisotropic distribution function (Ta = 10T/9, 
with T the standard two-body relaxation time) . Our value 
of Aa = 0.1 was fitted to averaged iV-body simulations of 
N < 2000. Finally Cb is the numerical parameter scaling 
the energy generation by three-body binaries, whose stan- 
dard value (Goodman 1987) is Cb = 90, which we take here. 
Goodman & Hut (1993) have recently argued that the for- 
mation rate of three-body binaries had been overestimated 
in previous work, so Cb = 75 would ensue for their new re- 
sults. But the still adopted value of Cb = 90 is within the 



Figure lb. Same quantities as in Fig. la, but per iV-body time 
unit per particle. 

expected error range, for their derivations as well as the ac- 
curacy is concerned with which we can fix Cb by comparison 
with the TV-body simulation. The reader more interested in 
the details of how to define parameters like A, Aa, and Cb 
and how to select their values is referred to GS and the other 
cited papers above. 

Fig. 2a compares the evolution of Lagrangian radii 
containing 1-90 % of the total mass between AG and our 
TV-body model. Whereas the pre-collapse evolution agrees 
extremely well, and also the slope of the post-collapse ex- 
pansion agrees fairly well, there is a difference in the col- 
lapse time itself, which is t — 2080 for the gaseous model 
and t = 2380 in the TV-body system. So is there something 
wrong? One has to recall that the conductivity parameter 
A in the gaseous model was adjusted to match simultane- 
ously the collapse time and slope of the Lagrangian radii 
evolution in pre-collapse for averaged TV-body simulations 
(N = 250, 500, 1000) and isotropic Fokker-Planck models 
(GS). In these averaged models there was always a consider- 
able spread in collapse times of individual TV-body models, 
as can be seen from Table 1. First we show in this table 
results obtained from various independent TV-body calcula- 
tions (TV < 2000), as e.g. the variance of Lagrangian radii for 
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the initial model or, minimum and maximum collapse time 
(£ m in and imax) as measured in individual simulations, the 
average collapse time t cc for all simulations of a given parti- 
cle number, the variance of the core collapse time a cc and its 
relation to the collapse time itself. For more details about 
the averaging process and the individual simulations see GHI 
and GHII, from whose models these data were taken. For all 
quantities except i m i n and t ma x there is an expected scaling 
with N denoted in Table 1; we found it in general agreement 
with the measured values, although there is some tendency 
for mismatch with such scaling for low N (N < 500). In 
particular the scaling law for the collapse time stems from 
the assumption that it is two-body relaxation (by taking 
a Coulomb logarithm factor of 7 = 0.11, see GHI, GHII) 
which dominates the evolution to core collapse. 

We have used these scaling laws to extrapolate quanti- 
ties for N = 10000 and N = 20000, which cannot be mea- 
sured yet, because there are not enough individual A"-body 
simulations (values for N = 1000 were used to extrapolate 
for N = 10000, and those of N = 2000 to extrapolate to 
N = 20000). We find for 10000 particles an extrapolated 
average collapse time of t cc = 2215. Note that even this 
value is not very accurate, because the measured average 
collapse times do not exactly obey the assumed relation. 
But if we take this value for the time being as a fiducial 
value it follows, that our individual iV-body simulation has 
a collapse time 165 time units (1.14 <r) larger than the aver- 
age, and gaseous and Fokker-Planck models have a collapse 
time 135 time units (0.93 a) smaller than the expected av- 
erage collapse time. 

We stress these numbers here in detail, because recently 
it has been claimed that the collapse time in anisotropic sys- 
tems could be somewhat larger as in the isotropic case, as 
suggested by new 2D Fokker-Planck models of Takahashi 
(1995, 1996), and previously by higher order moment mod- 
els of Louis (1990). However, at this present particle number 
of N = 10000 it seems very difficult to distinguish such a 
variation from the intrinsic uncertainties in determining the 
core collapse time. There are several factors which influ- 
ence measured collapse times, as e.g. in iV-body models 
the number of independent simulations, and in gaseous and 
Fokker-Planck models variation of quantities like tto and Cb 
(starting time and strength of binary energy generation, see 
GS for their definition, and their Figs. 4 and 5 for their influ- 
ence on collapse times) and the factor 7 used in the Coulomb 
logarithm could always cause a few per cent change in the 
collapse time. For example, 7 = 0.11 has been adopted as 
the best value from GHI and GHII, who considered only 
A" < 2000; a smaller 7 could always cause a larger collapse 
time. Hence we conclude here that from the present data 
there is some suspicion that the collapse time in iV-body 
simulations might be longer than expected from anisotropic 
gaseous and isotropic Fokker-Planck models, however, we 
consider the difference of just one a as not yet very con- 
vincing. Note that the free parameter in the conductivity of 
the anisotropic gaseous model was scaled to the same col- 
lapse time as the isotropic Fokker-Planck model (GS). Thus 
if one came to the conclusion that anisotropic 2D Fokker- 
Planck models and naturally anisotropic iV-body systems 
have a slightly longer collapse time, this just would mean 
that the anisotropic gaseous model had to be rescaled in 
its conductivity parameter A, now with respect to the 2D 



Fokker-Planck model. To judge about these questions has 
to be postponed to the future, when iV-body simulations of 
large particle numbers (N > 10000) become available in a 
larger number. 

In Fig. 2a one can see at a time much later than core 
bounce (t « 4000) that the gaseous model becomes unsta- 
ble to gravothermal oscillations at very late times (compare 
for gravothermal oscillations Bettwieser & Sugimoto 1984, 
Goodman 1987, Breeden, Cohn & Hut 1994). Since the in- 
stability at N = 10000 is not very strong (particle number 
not much above the critical number 7000 for gravothermal 
oscillations in the equal mass system) the oscillations start 
so late, at a time, which our iV-body simulation did not 
yet reach. It exhibits, however, some kind of oscillations 
directly after core bounce. But, as will become clear dur- 
ing the discussion of the binary effects, we cannot establish 
firmly that these oscillations are gravothermal, in contrast 
to just directly driven by binary energy. Gravothermal os- 
cillations (or expansion) should continue even after binary 
energy production has ceased. 

To relate our single N = 10000-body model with the 
lower N cases of GHI and GHII we have scaled its time 
coordinate by C(No)/C(N), where C(x) := x/log^x) (7 = 
0.11) for N = 10 4 and N = 1000, and plotted the data 
for the Lagrangian radii together with those of an averaged 
A" = 1000 body calculation, kindly made available by M. 
Giersz, in Fig. 2b. Such scaling would occur if standard two- 
body relaxation is the dominant force of evolution, which in 
fact can be clearly seen in the plot. The separation of the 
curves for different N at core bounce is due to the three- 
body activity which sets in at different times and scales with 
a different power of N. 

In order to compare the post-collapse evolution with the 
gaseous model it is better to scale the times of the gaseous 
model to the exact collapse time (t = 2380) of the iV-body 
simulation. The result on an enlarged timescale can be seen 
on Fig. 3. To get Fig. 3 we have multiplied the time and 
all radii with a free factor, so as to match the collapse time 
and the vertical position of the gaseous model and iV-body 
curves. Consistently with the previous results of GHI, GHII, 
and GS we find that core bounce in the iV-body model is 
deeper in the inner shells and that the intermediate zones 
of the gaseous model expand faster, whereas the outermost 
shells lag behind as compared to the TV-body model. We 
think this is due to a non-local energy transport by high- 
velocity stars on radial orbits, which are reaction products 
from superelastic binary encounters in the core, as elabo- 
rated also in the previously cited papers. In the gaseous 
model energy can only be transported locally by heat con- 
duction with the temperature gradient; thus the energy cre- 
ated by binaries, which is of the correct magnitude (see 
below) causes more expansion in the adjacent, intermedi- 
ate shells relative to the iV-body model. The importance 
of high-velocity stars carrying energy from superelastic en- 
counters quickly out of the core was first realized by Aarseth 
(1974) and further elaborated on by Makino & Sugimoto 
(1987), who denoted them as suprathermal particles. 

Figs. 4a-d show the radial and tangential velocity dis- 
persions, each for the entire evolutionary time (4a, 4c) and 
for an enlarged time segment around core bounce (4b, 4d). 
In most cases the general agreement between the gaseous 
model and our single A"-body simulation is striking. How- 
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Table 1. Radial variance of initial models an, mini- 
mum and maximum core collapse times t m j n and t max , 
average core collapse time t cc , variance of core collapse 
time (Tcc, and the ratio cr cc /t cc as a function of particle 
number N. Values with a *-symbol are estimated from 
the scaling laws indicated in the last lines. All data for 
N < 2000 kindly supplied by M. Giersz (priv. comm.) 
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Figure 3. Lagrangian radii as in Fig. 2a, but now an enlarged 
time segmented near core bounce. Gaseous model rescaled ac- 
cording to text in order to compare the post-collapse slopes of 
both models. The times of prominent hard binary events are 
marked, as there are binary escapers (dots) and strog hardening 
events of binaries by three-body encounters (stars). 



Figure 2a. Lagrangian radii containing the indicated fraction 
of total mass as a function of time for the 10 4 -body simulation 
(fluctuating curves) and the standard anisotropic gaseous model 
(smooth curves). 



Figure 2b. Lagrangian radii as in Fig. 2a, but here our 10 4 -body 
simulation (time rescaled according to main text) compared to an 
averaged N = 1000 simulation of GHI and GHII. 

ever, there are some features near and after core bounce 
worth discussing in more detail. First, the effect that inter- 
mediate shells (e.g. 10-20 % of total mass) expand a little 
faster than the iV-body system, as was detected for the La- 
grangian radii, corresponds to a smaller "temperature" (ve- 
locity dispersion) there. On the other hand, further outside, 



Figure 4a. Radial ID velocity dispersion averaged between the 
indicated Lagrangian radii as a function of time. Comparison 
between Af-body (fluctuating curves) and gaseous model (smooth 
curves). 

where the Af-body system expands more rapidly, it is some- 
what cooler than the AG model. In general one could say 
that the large scale Figs. 4a and 4c support our claim of the 
generally fair agreement between A^-body and AG models, 
but Figs. 4b and 4d presenting a zoom in on core bounce and 
post-collapse give a clearer picture of what are the remaining 
differences. How about the anisotropy A = 2— 2of /<rj? itself? 
In principle it is redundant, all information is already con- 
tained in the previous plots. However, the analysis of A en- 
hances the remaining differences between the models, so we 
will discuss it here as well with the help of Fig. 4e. For the 
innermost Lagrangian radii (say up to 30 %) it is not useful 
to discuss A, because it is fluctuating strongly; comparisons 
with the AG models have to rely on the velocity dispersions 
as above. Note, however, the good agreement between both 
models for the anisotropy averaged between 30 and 40 % 
of total mass (lowermost curves of Fig. 4e), especially in 
the reduction of anisotropy after core bounce. Surprisingly, 
there is not such agreement for the larger radii. For 40-50 % 
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Figure 4b. As Fig. 4a, but for an enlarged time section near 
core bounce. 



Figure 4e. As Fig. 4a, but for the anisotropy defined as 



Figure 4c. As Fig. 4a, but for the tangential ID velocity dis- 
persion. 



Figure 4d. 

persion. 



As Fig. 4b, but for the tangential ID velocity dis- 



the decrease of A in post-collapse is only visible in the AG 
model, the anisotropy in the Af-body model increases fur- 
ther. More outside, the disagreement between both mod- 
els starts even earlier. This is consistent with the idea of 
suprathermal particles, moving on radially elongated orbits 
towards the halo, not being accounted for in the gaseous 



model. As for the outermost shells of the system one has 
also to take into account a difference in the boundary con- 
ditions (the A-body model is open and may expand freely 
into space, the AG model is confined by an adiabatic wall at 
100 % of total mass, so the 100 % Lagrangian radius cannot 
expand, which it does in the A-body model). 

Figs. 5a-c depict the time development of the central 
potential, Cohn's x value (x = $ c /3cr;?, where a c is the core 
velocity dispersion), and the number of core particles. All of 
the data in Figs. 5 are scaled in time as in Fig. 2b, in order 
to compare with the data of GHI and GHII. In some cases 
the figure caption notes that the data were smoothed, which 
means the high frequency fluctuations have been cut out 
by the "Numerical Recipes" routine SMOOFT, which uses 
Fourier transforms to do that. The window for smoothing 
was about 30 data points, which are 45 time units. In Figs. 
5a and c once more a very good agreement of a single N = 
10 4 model with the average N = 1000 data is visible for the 
pre-collapse phase. During core bounce the single N = 10 4 
simulation exhibits strong fluctuations (very deep central 
potential, correlated with an extremely small core particle 
number), which presumably would be smoothed out if we 
would be able to average the results for different independent 
models, as in the N = 1000 case. 

It is very instructive to look at a graph just showing as 
tiny dots values of the potential at the individual particle 
locations as a function of radius for different times in pre- 
and post-collapse (Fig. 6). This gives a good feeling how 
smooth the particle distribution already is for a total parti- 
cle number of 10000, and moreover, it unambigously clarifys 
that core bounce, the time with the deepest central poten- 
tial, occurred at t = 2380 (labelled 886 in the plot), which 
can also be seen in Fig. 5a (but note the rescaled time) . 



3.2 Escapers and binaries 

At the present endpoint of our simulation a total of 360 par- 
ticles had escaped (they were removed from the calculation 
when speeding with more than the escape velocity from the 
entire cluster and being at a radius larger than ten times 
the initial scaling radius of Plummer's model). Fig. 7a,b 
show the cumulative particle and energy loss. Similarly to 
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Figure 5a. Central Potential as a function of scaled TV-body 
time, see text. Compared with a smoothed curve obtained from 
average TV = 1000 data. 



Figure 5c. Core particle number as a function of scaled TV-body 
time. Compared with data of average lower TV simulations. 



Figure 5b. Cohn's x as a function of scaled TV-body time. 
Smoothed in time, compared with smoothed data from average 
TV = 1000 simulation. 

GHI we tried to approximate the TV-body results by rescal- 
ing Henon's (1965) escape rate (isotropic Plummer sphere) 
multiplicatively for the effects of evolution and anisotropy 
(see Figs. 9 and 11 of GHI). The effect of evolution was 
determined by directly calculating the escape rate of an 
isotropic Fokker-Planck model with TV = 10000 particles and 
standard parameters (7 = 0.11 for the Coulomb logarithm, 
C\ — 90 for the energy generation), see for reference Eq. 2 
of GHI. As a model for the effect of anisotropy Dejonghe's 
(1987) anisotropic Plummer models were taken, which have 
a dimensionless parameter q. q can be determined by tak- 
ing the anisotropy A at a certain radius r of the TV-body 
model by q — A/(l-\- 1/V 2 ). GHI took the anisotropy at 
the Lagrangian radius containing 75 % of total mass, here 
we try to bracket the TV-body results by computing q and 
subsequently scaled escape rates for both the 75 % and 90 % 
Lagrangian radius (see plots). Since we were able to distin- 
guish every single escaper in our TV-body calculation, it was 
possible to remove those escapers by number and energy, 
which are due to three-body encounters (scatterings with 
hard binaries, in the following denoted as "hard escapers"), 
defined by the condition that their energy of escape is larger 
than 3 kT. Such reduced particle and energy loss is labelled 



Figure 6. Individual particle's potentials for different times in 
pre- and post-collapse as indicated. 

by N-BIN in the plots. 

From the figures we conclude the following: the effect 
of hard escapers in the TV-body rates of escape is small for 
the particle number and in pre-collapse, but dominant for 
the energy of escapers in core bounce and beyond. A fair 
approximation of the particle and energy loss in the TV-body 
system is reached by rescaling Henon's escape rate for the 
effect of evolution and anisotropy, provided the energy of 
the hard escapers is subtracted first from the TV-body data. 
The effect of evolution visible in the curves labelled FP-I 
is much less significant than that of the anisotropy. How- 
ever, since Dejonghe's models are only a very approximate 
model it is not a priori clear, how to fit the best q parame- 
ter to an individual TV-body system. To illustrate this un- 
certainty we have used two different Lagrangian radii (con- 
taining 75 % and 90 % of total mass, with correspondingly 
different anisotropies and thus varying q parameter for the 
same time). The escape rate of particles seems to be best 
approximated by using a g-value derived from the 90 % La- 
grangian radius (Fig. 7a), whereas the energy loss rate is 
bracketed by the two choices (Fig. 7b). It follows that it 
is not possible to find a simultaneous best fit to both data; 
we think that this is due to the fact that Dejonghe's models 
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are not an accurate representation of the real iV-body distri- 
bution. Note that we had to rescale here once more the N- 
body time to the same collapse time as in the Fokkcr-Planck 
model (t — 2080, see above Fig. 2a and its discussion) used 
for comparison, in order to get a useful comparison of both 
data. 

We started initially without any binaries. At t = 1200 
a transient mildly hard binary formed, whose effect can be 
seen in a slight dip of the iV-body Lagrangian radii against 
the gaseous model data in Fig. 2a. Such very early transient 
binary activity may be a reason, why our model system has 
a collapse time, which is at more than one a above the ex- 
pected average value. The really hard binaries then formed 
shortly before core bounce. This can be seen on Figs. 8a, 
b, where the high energy (E >• lkT) escapers occur just 
at core bounce, which is the time of formation of hard bi- 
naries, as is shown in Fig. 8c. Fig. 8b is an enlargement 
from the data of Fig. 8a; we see that the binary activ- 
ity reduces to a significantly lower level again at t — 2600, 
since we interpret any escaper with more than lkT as due 
to an interaction with a hard binary. It is interesting to 
note, that near the times t = 2400, 2500, and 2600 there 
a crowding of high energy escapers in Fig. 8a. We relate 
this to the spikes in the velocity dispersions of the outer 
Lagrangian shells occurring in Figs. 4 just near the same 
time, only for the radial velocity dispersion, not for the tan- 
gential one. This supports the idea, that the high-energy 
reaction products of the strong three-body encounters move 
radially outwards and generate the spikes in oy. Moreover, 
these three times can be identified as periods of strong bi- 
nary activity in Fig. 8c (binaries exceeding a binding energy 
of 200 in iV-body units). If one then looks at Fig. 3 these 
times approximately correspond to expansion phases of the 
Lagrangian radii; after the binary activity has ceased, the 
expansion does not continue, on the contrary, the system 
starts to contract again. Note, that after t — 2600 there is 
also no strong binary activity, because the hardest binary, 
although it is still contained in the system (not yet escaped), 
moves on an elongated orbit after a recoil. Hence it has a 
small probability for three-body encounters. From all this 
evidence we conclude that the post-collapse oscillations of 
the iV-body system are binary driven and not gravother- 
mal. If they were gravothermal, expansion should continue 
after the binary energy generation has ceased. Moreover, as 
already Bettwieser & Sugimoto (1984) stated there should 
be a radial temperature (i.e. velocity dispersion) inversion 
in the case of a gravothermal expansion, which we could 
not observe with a sufficient statistical significance in our 
iV-body data. But on the other hand, there is a steady ex- 
pansion consistent with the gaseous model in post-collapse 
superposed on the binary driven oscillations. So we would 
not exclude that on a secular timescale our TV-body system 
will undergo a gravothermal expansion, presumably with a 
temperature inversion too small to be extracted from the 
noisy data. 

In Fig. 8c the fate of the hard binaries in the system can 
be followed individually, e.g. one can estimate by eye that 
the rate of change of the binding energy of a hard binary 
is related to its binding energy x; in an encounter binaries 
with large x suffer a change in binding energy of the same 
order as x. This is consistent with the early experimental 
(Hills 1975) and theoretical (Heggie 1975) predictions, that 



Figure 7a. Cumulative particle loss due to escapers as a function 
of scaled TV-body time (scaled to the collapse time of an isotropic 
Fokkcr-Planck model). Curves labelled N and N-BIN denote the 
TV-body data with and without hard escapers (see main text); the 
data labelled FP-I show the increase of Henon's (1975) standard 
escape rate (at t = 0) due to the evolutionary effect; the curves la- 
belled A-90 and A- 75 show the multiplicative combination of the 
evolutionary effect and another correction for anisotropy evalu- 
ated by using Dcjonghe's (1987) anisotropic Plummer model with 
a g-value derived from the anisotropy at the 75 % or 90 % La- 
grangian radius, respectively. 



Figure 7b. As Fig. 7a, but for the energy of the escapers. 

on average the energy change Ax during a superelastic scat- 
tering should be Aa; = ijx (n = 0.4). GHII argue that an 
experimental value of rj, found by their statistics of a large 
number of TV-body simulations for TV < 2000, is about 0.2, 
in better agreement with more recent predictions by Spitzer 
(1987) and scattering experiments by Heggie & Hut (1993). 
Since, however, in contrast to the averaged TV-body simula- 
tions of GHI and GHII, we only have a rather small number 
of strong encounters, i.e. a poor statistics, we will not elab- 
orate further on this point here. We think our results are 
at least consistent with a value of rj between 0.2 and 0.4; 
a change with r\ = 0.4 is indicated in the figure as a thick 
vertical bar for three energies x — 20, 100, 300. One may 
wonder why in one case at about t — 2570 the binding en- 
ergy of a binary decreases so dramatically. This is due to the 
formation of a hierarchical triple, which was unfortunately 
not included in our binary list. Later, its binding energy 
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Figure 8a. Time at which stars escape as a function of their en- 
ergy in kT, where lkT = with the ID core velocity dispersion 

CTc- 



Figure 8c. Individual fate of binaries close to core bounce; bind- 
ing energy as a function of time. At some times labels are plotted, 
which are the names of the stars in the binary and can be used to 
identify them in the following plot. Three thick vertical line seg- 
ments show the magnitude of an increase in binding by energy by 
40% from energy 20, 100, and 300, respectively, to be compared 
with the actual energy changes. Binary escapers and exchange 
reactions of the binaries are indicated. 



Figure 8b. As Fig. 8a, but with an enlarged time axis to see 
the decrease in binary activity after core bounce for t > 2600. 

occurs again in the most strongly bound binary, after disso- 
lution of the triple. Fig. 8d gives an optical impression of 
the orbital motion of the centres of masses of hard binaries, 
their recoils by superelastic scatterings to very large radii, 
and subsequent shrinking again by dynamical friction to the 
core. Especially the motion of the hardest binary consisting 
of the particles # 4805 and 5745 (see Fig. 8c, its binding 
energy amounts to several hundred kT) reaches far into the 
halo after an energectic recoil. Thus one can understand 
that exotic products of binary evolution (like soft X-ray bi- 
naries or blue stragglers, see for the first e.g. Rappaport et 
al. 1994) may be found in the outskirts of a cluster. 

As a final remark we want to demonstrate with Fig. 8e 
how the cumulative binary energy generation in the A^-body 
system compares to the steady energy input modelling the 
binary heating in the gaseous model. One can see that the 
energy liberated in the gaseous model is of the same order as 
the binding energy of the hard binaries still in the Af-body 
system, and also of the same order as the binding energy 
carried away by escaping binaries. 

3.3 Movement of the density centre 

A very interesting and not yet completely understood fea- 
ture of a real Af-body system is the oscillations of its density 



Figure 8d. Radial position of the hard binaries of the previous 
figure over the same time span. Note that the core radius is 
much smaller than one during this entire period. In some cases 
by comparison with Fig. 8c the event can be identified, which 
kicked a binary out of the core or even out of the cluster (lines 
escaping to the top belonging to escaping binaries). 

centre and also other quantities, like the Lagrangian radii. 
To define the density centre we first compute for each parti- 
cle i a neighbour density pi, according to Casertano & Hut 
(1985) by the mass inside r& (excluding the test particle it- 
self), where r& is the distance of the sixth nearest neighbour. 
Then by 

(0 (») 

we define a density centre for a special subset of particles 
(r d , Vi are the position vectors of the density centre and the 
particle i, resp.). For the summation only the inner parts 
of the system should be included, in order to exclude strong 
biases by asphericities in the outermost parts of the system 
(large weight due to large | | ) . On the other hand selecting 
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Figure 8e. Energy balance as a function of time near core 
bounce. Labelled are the curves for the binding energy of the bi- 
naries still inside the cluster (Ebin): the cumulative translational 
energy of single cscapcrs (E esc ), the cumulative binding energy 
carried away by escaping binaries (E e scbin)> and the cumulative 
energy due to three-body binaries produced in a standard AG 
model (-Ebingas)- The core collapse time of the gaseous model 
was scaled to the actual collapse time of the N-bodj model, see 
text. 

the core particles or some fraction of them (based on the 
previous core radius, which was also determined in a sim- 
ilar way as the density centre) proved to produce extreme 
fluctuations of the density centre in the late phases of core 
collapse. So we chose a procedure suggested by McMillan, 
Hut & Makino (1990), to take all particles whose pi is some 
fraction (actually 1/20) of the maximum pi. This yielded 
the data for the density centre of Figs. 9a-c. In recent 
years there were detailed examinations of core wandering 
and oscillations by Heggie, Inagaki & McMillan (1994) and 
Sweatman (1993). Our data are different from their data in 
several respects. Firstly, we cover a much larger time, how- 
ever, do not resolve the small timescales as they did. We 
cover approximately 1200 N-body time units with a con- 
stant time interval of roughly half an initial mass crossing 
time (At = 1.4). In principle our dataset is even larger, but 
unfortunately we did not always keep a constant time inter- 
val for the data collection, which turns out to be essential for 
the subsequent analysis. Secondly, we are not able to resolve 
the noise at high frequencies, as was the aim of Heggie, Ina- 
gaki & McMillan (1994). In fact, our data become significant 
just at frequencies smaller than 0.1 (see their Fig. 7), where 
their statistics becomes poor. Also Sweatman (1993), who 
did a very thorough examination, which we closely follow 
here in our data processing, was limited to periods smaller 
than 10 (frequencies larger than 0.1), due to the duration of 
his simulation. 

We first present the original data of the density centre in 
three different scales in Figs. 9a-c, together with the centre 
of mass data for the entire system. As noted already by the 
cited authors and also by Makino & Sugimoto (1987) there 
are variations on different timescales, beginning from very 
abrupt changes on times comparable to the dynamic time 
up to larger changes on timescales of the order of several 
crossing times. We will now inspect the time variability of 
the data in more detail. 



In Figs. 9b, c a third (thick) curve is plotted, which 
is the density centre data r s d smoothed by the SMOOFT 
procedure of "Numerical Recipes" with a window of 50 data 
points (we varied the window to 30 and 100 data points with- 
out any significant change in the following results) . The rea- 
son is that for a proper autocorrelation and frequency anal- 
ysis we need data whose average is zero. Sweatman (1993) 
took the centre of mass as a reference point, but as one can 
see in Fig. 9a for a long time simulation, where many, partly 
high energy, escapers are removed and the outermost parts 
of the system are not nicely symmetric anymore, there is a 
systematic shift between the centre of mass and the density 
centre. Therefore we need the smoothed data as reference, 
such that we can compute a quantity SR — r s d — r<j, whose 
time average is zero. Therefore we cannot detect any periods 
in the data which are larger than the smoothing window, so 
we checked that our results are independent of the size of it. 
With SR = |5R| we performed an autocorrelation analysis 
identical to that described by Sweatman (1993) in his chap- 
ter 3.3. The striking result is in Fig. 9d, which exhibits a 
dominant period of the order of 40 time units, which is about 
14 half-mass crossing times. To check this result indepen- 
dently we took the data for SR and performed numerically 
a discrete Fourier transform (by use of the Four-function of 
Mathematica). The resulting power spectrum as a function 
of period is displayed in Fig. 9e - again there is a clear max- 
imum near periods of the order of 50, one could distinguish 
two subpeaks at about 40 and 60. 



4 CONCLUSIONS AND DISCUSSION 

We have performed a very long (in terms of computer time as 
well as of physical time) simulation of a collisional N = 10 4 
particle system well into the post core collapse evolutionary 
phase by using the standard NBODY5 integrator (Aarseth 
1985). As initial model a Plummer sphere with no bina- 
ries and equal masses was selected. In contrast to previous 
averaged AT-body simulations (Giersz & Heggie 1994a, b, 
GHI, GHII; Giersz & Spurzem 1994, GS) this single model 
already compares well with expectations of gas dynamical 
models based on the Fokker-Planck approximation. In par- 
ticular the evolution until core collapse is entirely domi- 
nated by standard two-body relaxation, in particular La- 
grangian radii, velocity dispersions and anisotropy compare 
fairly well with the other, statistical models, although there 
are some disagreements, especially in post-collapse and for 
the outer halo regions, which we think are due to the poor 
representation of suprathermal particles, originating from 
superelastic scatterings with hard binaries, in the gaseous 
model. Anisotropy is also the dominant factor determining 
the time evolution of the escape rate, which in deep core 
collapse is much larger than Henon's (1965) standard esti- 
mate. Within some uncertainty the effect of the anisotropy 
on the escape rate can be modelled by using one of De- 
jonghe's (1987) anisotropic Plummer models, provided the 
energy of escapers from hard three body encounters is not 
taken into account. 

Thereafter the core bounces due to the energy gener- 
ation by superelastic scatterings between single stars and 
binaries. The energy generation is in fair accord with statis- 
tical predictions by Heggie (1975) and Hut (1993). From our 
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Figure 9a. Modulus of the vector of the density centre and the 
centre of mass of the entire cluster as a function of time. 



Figure 9d. Autocorrelation Function explained in the text as a 
function of displacement (in time) A (in N-body time units). 



Figure 9b. As Figure 9a, with enlarged time axis, and added 
the thick curve for the smoothed density centre variations, which 
serves as a reference for the autocorrelation analysis, see text. 



Figure 9c. As Figure 9a, but even more zoomed in on the time 
axis, to see the limiting resolution of our data at small time scales. 

simulation we see that individual hard binaries are some- 
times ejected by recoil effects very far from the core after 
encounters. Since in our simulation we are able to follow 
the individual fate of each binary, we can correlate strong 
binary activity with expansion phases of the whole system. 



Figure 9e. Power spectrum of the differences between the den- 
sity centre data in Figs. 9a-c and their smoothed counterpart (see 
Figs. 9bc), as a function of period in TV-body time units. 

Since such expansion does not continue after the binary ac- 
tivity ceases, the post-collapse oscillations present in our 
model are not gravothermal in nature. Although theoreti- 
cally 10000 particles should be enough to observe gravother- 
mal oscillations (but they set in long time after core bounce 
in the gaseous model), the fluctuations induced by individ- 
ual strong three-body encounters suppress them to a non- 
observable level. Nevertheless the A^-body expands on a 
secular timescale as prescribed by the gaseous model, so 
it cannot be excluded that there is a gravothermal expan- 
sion with superimposed binary-driven oscillations. In recent 
simulations (Makino 1996) with even larger particle number 
(32k particles) the onset of gravothermal oscillations includ- 
ing a temperature inversion has been observed, which would 
mean that for increasing N the gravothermal evolution as 
predicted by gaseous and Fokker-Planck models will domi- 
nate the evolution as compared to the stochastic variations 
directly induced by binary activity. For a particular case, 
namely the fluctuations of core collapse time in an individ- 
ual A/-body simulation, we could show that it will decrease 
relative to the collapse time itself with larger N. 

Most features of our single N = 10 4 model are consis- 
tent with expectations from statistical (gaseous or Fokkcr- 
Planck models), although our collapse time is 1.14 a larger 
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than an extrapolated average value for iV-body simulations, 
and the latter is another 0.93 a larger than the collapse time 
of anisotropic gaseous and isotropic Fokker-Planck models. 
Since this is rather large, one may speculate whether it re- 
flects a recent result by Takahashi (1995, 1996), who finds 
that in his new 2D FP models collapse times are some 10 
% larger than expected from the previous isotropic models. 
However, our particle number is still not large enough to 
discriminate with a sufficient statistical significance between 
such small differences in the collapse time. 

We have also examined the wandering of the density 
centre motions. By autocorrelation analysis and fourier 
transformation we find a dominant period of the order of 40- 
60 time units, which is about 13-20 initial half-mass crossing 
times. Since our database is much longer in time, these os- 
cillations could not have been found in previous studies like 
Heggie, Inagaki & McMillan (1994) and Sweatman (1993). 
The existence of long-time scale oscillations like those we 
found here has, however, been noted already by Makino & 
Sugimoto (1987). 

From data of our simulation which were not shown in 
the plots here we conclude that the velocity distribution in 
the core changes significantly between the pre- and post- 
collapse state of the cluster. Thus for example the ratio of 
the fourth order moment of the velocity distribution related 
to the second order moments is larger after core bounce than 
before, consistent with the idea that a high velocity tail of 
stars is generated by superelastic scatterings with the bi- 
naries. The presentation of these data as well as possible 
relations to observable line profiles in star clusters will be 
subject of future work. In parallel to this it could be useful 
to extend the gaseous models to higher order moments of 
the velocity distribution in order to achieve a better repre- 
sentation of the suprathermal particles in the post-collapse 
phase. One might argue that 2D Fokker-Planck simulations, 
if they become more common in the near future, already 
deliver a complete and unrestricted 2D distribution func- 
tion. However, for the grid-based solutions of the orbit- 
averaged Fokker-Planck equation still an isotropized distri- 
bution function is used to calculated the diffusion coeffi- 
cients. This is in contrast to the gaseous models, which use 
a consistent anisotropic background for the Fokker-Planck 
collisional terms. If the Fokker-Planck equation is solved by 
using a variational principle (Inagaki & Lynden-Bell 1990, 
Takahashi & Inagaki 1992) it is self-consistent (contrary to 
the discretization method) , but there occur other uncertain- 
ties, as e.g. the proper choice of trial functions. We con- 
clude that none of the presently used methods to simulate 
the evolution of large N star clusters should be discarded; 
it will remain very important in the future to produce a 
larger data base, especially for the case of iV-body simula- 
tions with large N and to thoroughly study the strong and 
weak features of each method. 

To improve the statistics at particle numbers of N > 
10000 more independent computations are useful and in 
progress, especially by using the recent fast versions of the 
HARP special purpose hardware (Taiji 1996); to see whether 
a single star cluster behaves in accord with the statistical 
models and, for example, shows gravothermal oscillations, a 



still larger N is required, as in the case of the recent 32000 
particle simulation (Makino 1996). 
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